Wczytanie danych

# wczytanie pakietu
library("terra")

W pierwszym kroku musimy wylistować dane (rastry), które zamierzamy wczytać. W tym celu możemy wykorzystać funkcję list.files(), która jako argument przyjmuje ścieżkę do folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy wczytać (pattern = "\\.TIF$") oraz zwrócić pełne ścieżki do plików (full.names = TRUE).

# listowanie plikóW z katalogu
files = list.files("dane/landsat", pattern = "\\.TIF$", full.names = TRUE)
files
## [1] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF"
## [2] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF"
## [3] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF"
## [4] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4.TIF"
## [5] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5.TIF"
## [6] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6.TIF"
## [7] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7.TIF"

Kiedy utworzyliśmy już listę plikóW, możemy je wczytać przy pomocy funkcji rast() z pakietu {terra} i następnie wyświetlić metadane.

# wczytanie danych rastrowych
landsat = rast(files)
landsat # odwołanie się do obiektu wyświetla metadane
## class       : SpatRaster 
## dimensions  : 1853, 2846, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 602535, 687915, 5742525, 5798115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
## sources     : LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF  
##               LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF  
##               LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF  
##               ... and 4 more source(s)
## names       : LC08_~SR_B1, LC08_~SR_B2, LC08_~SR_B3, LC08_~SR_B4, LC08_~SR_B5, LC08_~SR_B6, ... 
## min values  :          21,          32,        1321,        1924,        6827,        6982, ... 
## max values  :       60441,       61519,       65451,       65192,       63873,       65454, ...

Możemy również skrócić lub zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.

names(landsat) # nazwy oryginalne
## [1] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1"
## [2] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2"
## [3] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3"
## [4] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4"
## [5] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5"
## [6] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6"
## [7] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7"
names(landsat) = paste0("B", 1:7) # skrócenie nazw
names(landsat) # nowe nazwy
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
# zamiana nazw
# names(landsat) = c("Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2")

Wczytanie danych wektorowych odbywa się w analogiczny sposób za pomocą funkcji vect().

# wczytanie danych wektorowych
poly = vect("dane/powiat_sremski.gpkg")
poly
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 1  (geometries, attributes)
##  extent      : 622510.6, 659710.4, 5754795, 5784772  (xmin, xmax, ymin, ymax)
##  source      : powiat_sremski.gpkg
##  coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
##  names       : TERYT
##  type        : <chr>
##  values      :  3026

Teraz możemy przygotować prostą wizualizację jednego kanału (np. ultra blue; B1) oraz poligonu.

# wizualizacja
plot(landsat[[1]])
plot(poly, add = TRUE)

Operacje na rastrach

Zasięg naszego obszaru analizy ograniczony jest do powiatu śremskiego, natomiast scena satelitarna ma o wiele większy zasięg. W takiej sytuacji możemy dociąć rastry, dzięki czemu przetwarzanie będzie szybsze, a finalna mapa zajmie mniej miejsca na dysku. Do docinania rastrów służy funkcja crop() i jako argumenty musimy podać raster oraz wektor.

landsat = crop(landsat, poly)
plot(landsat[[1]])
plot(poly, add = TRUE)

Obszar rastra zmniejszył się. Jednak możemy zauważyć, że poza poligonem wartości nie zostały usunięte. Wynika to z faktu, że obraz zawsze docinany jest do obwiedni (bounding box), a wartości poza obiektem/poligonem w rzeczywistości są maskowane (tj. są oznaczane jako brakujące wartości). Aby zamaskować piksele poza poligonem należy użyć funkcji mask().

landsat = mask(landsat, poly)
plot(landsat[[1]])
plot(poly, add = TRUE)

Tą operację można przeprowadzić również w jednej linii kodu używając argumentu mask = TRUE w funkcji crop().

crop(landsat, poly, mask = TRUE)

W następnym kroku możemy w prosty sposób sprawdzić statystyki opisowe naszego zbioru danych.

summary(landsat)
##        B1              B2              B3              B4       
##  Min.   : 6577   Min.   :   32   Min.   : 6844   Min.   : 7686  
##  1st Qu.: 8039   1st Qu.: 8220   1st Qu.: 9093   1st Qu.: 8649  
##  Median : 8425   Median : 8614   Median : 9609   Median : 9353  
##  Mean   : 8650   Mean   : 8921   Mean   : 9903   Mean   : 9961  
##  3rd Qu.: 9044   3rd Qu.: 9370   3rd Qu.:10396   3rd Qu.:10796  
##  Max.   :19736   Max.   :21314   Max.   :24431   Max.   :27674  
##  NA's   :48391   NA's   :48390   NA's   :48390   NA's   :48390  
##        B5              B6              B7       
##  Min.   : 7240   Min.   : 7378   Min.   : 7357  
##  1st Qu.:14669   1st Qu.:12482   1st Qu.:10002  
##  Median :17462   Median :14082   Median :11417  
##  Mean   :18004   Mean   :14673   Mean   :12519  
##  3rd Qu.:21254   3rd Qu.:16616   3rd Qu.:14170  
##  Max.   :31081   Max.   :50865   Max.   :60130  
##  NA's   :48390   NA's   :48390   NA's   :48390

Jak możemy zauważyć wartości odbicia spektralnego dla naszego zbioru danych są w zakresie od kilku do kilkunastu tysięcy dla każdego kanału. Z definicji odbicie spektralne jest w przedziale od 0 do 1. Takie dane musimy przeskalować za pomocą poniższego równania:

\[x = x \cdot 0.0000275 - 0.2\]

Nie ma potrzeby stosowania tego wzoru osobno dla każdego kanału w pętli, ponieważ domyślnie operacje matematyczne są stosowane dla wszystkich kanałów.

landsat = landsat * 2.75e-05 - 0.2
summary(landsat)
##        B1              B2              B3              B4       
##  Min.   :-0.02   Min.   :-0.20   Min.   :-0.01   Min.   :0.01   
##  1st Qu.: 0.02   1st Qu.: 0.03   1st Qu.: 0.05   1st Qu.:0.04   
##  Median : 0.03   Median : 0.04   Median : 0.06   Median :0.06   
##  Mean   : 0.04   Mean   : 0.05   Mean   : 0.07   Mean   :0.07   
##  3rd Qu.: 0.05   3rd Qu.: 0.06   3rd Qu.: 0.09   3rd Qu.:0.10   
##  Max.   : 0.34   Max.   : 0.39   Max.   : 0.47   Max.   :0.56   
##  NA's   :48391   NA's   :48390   NA's   :48390   NA's   :48390  
##        B5              B6              B7       
##  Min.   :0.00    Min.   :0.00    Min.   :0.00   
##  1st Qu.:0.20    1st Qu.:0.14    1st Qu.:0.08   
##  Median :0.28    Median :0.19    Median :0.11   
##  Mean   :0.30    Mean   :0.20    Mean   :0.14   
##  3rd Qu.:0.38    3rd Qu.:0.26    3rd Qu.:0.19   
##  Max.   :0.65    Max.   :1.20    Max.   :1.45   
##  NA's   :48390   NA's   :48390   NA's   :48390

Nadal możemy zauważyć, że pewne wartości przekraczają nasz zakres od 0 do 1. Są to wartości odstające, które zazwyczaj związane są z błędnym pomiarem lub nadmierną saturacją. Można ten problem rozwiązać na dwa sposoby:

  1. Zastąpić te wartości brakiem danych (NA).
  2. Dociąć do minimalnej i maksymalnej wartości.

Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.

# sposób nr 1
landsat[landsat < 0] = NA
landsat[landsat > 1] = NA
# sposób nr 2
landsat[landsat < 0] = 0
landsat[landsat > 1] = 1

Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym przypadku zamiast funkcji plot() należy użyć funkcji plotRGB() oraz zdefiniować kolejność kanałóW czerwonego, zielonego oraz niebieskiego. Oprócz tego, należy podać maksymalną wartość odbicia dla kanałów (w naszym przypadku scale = 1). Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować rozciągnięcie histogramów używając argumentu stretch = "lin" lub stretch = "hist".

plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1, stretch = "lin")

Klasteryzacja

# wczytanie pakietu do klasteryzacji
library("cluster")

Modele do klasteryzacji (oraz klasyfikacji nadzorowanej) wymagają macierzy lub ramki danych. Zatem dane do treningu musimy przygotować w odpowiedni sposób. Raster można sprowadzić do macierzy przy użyciu funkcji values().

mat = values(landsat)
nrow(mat) # wyświetla liczbę wierszy
## [1] 1238760

Za pomocą funkcji View() możemy sprawdzić jak wygląda nasza macierz. Jak łatwo można zauważyć, mnóstwo jej wartości oznaczonych jest jako brak danych (głównie są to wartości poza obszarem analizy). Zazwyczaj modele nie obsługują NA, więc musimy je usunąć. Służy do tego dedykowana funkcja na.omit().

mat_omit = na.omit(mat)
nrow(mat_omit)
## [1] 640802

Teraz przejdziemy do najważniejszego etapu analizy, czyli do wytrenowania modelu. Użyjemy najprostszego modelu grupowania k-średnich (k-means). Ten model wymaga jedynie, aby podać z góry liczbę klastrów/grup (argument centers). Jest to algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby analiza była powtarzalna musimy ustawić ziarno losowości – set.seed().

set.seed(1)
mdl = kmeans(mat_omit, centers = 5)

W wyniku powyższej operacji otrzymaliśmy m. in.:

  1. Obliczone średnie wartości klastrów dla poszczególnych kanałóW (mdl$centers).
  2. Wektor ze sklasyfikowanymi wartościami macierzy (mdl$cluster).

Wyświetlmy te obiekty:

mdl$centers
##           B1         B2         B3         B4        B5        B6         B7
## 1 0.02308493 0.02740529 0.04435163 0.04147431 0.1727330 0.1309803 0.07958614
## 2 0.02044635 0.02579641 0.05547638 0.03750665 0.4421821 0.1492823 0.07336897
## 3 0.08051378 0.09522015 0.13257100 0.16603768 0.2716112 0.3437177 0.32147881
## 4 0.03431384 0.04119050 0.07134604 0.06554308 0.3509450 0.2042120 0.12686408
## 5 0.04921918 0.05830302 0.08442430 0.09825275 0.2417769 0.2551757 0.19605182
head(mdl$cluster) # wyświetla pierwsze 6 elementów
## [1] 4 5 5 1 1 1

Oznacza to, że pierwszy wiersz (piksel) należy do grupy 3, drugi do grupy 2, trzeci do grupy 2, itd. Kolejnym etapem jest stworzenie mapy na podstawie otrzymanego wektora z klastrami.

Na początku musimy przygotować pusty wektor składający się z całkowitej liczby pikseli rastra. Można to sprawdzić za pomocą funkcji ncell(). W naszym przypadku jest to 1 238 760.

vec = rep(NA, ncell(landsat)) # przygotuj pusty wektor

Następnie musimy przypisać nasze klastry w wektorze w odpowiednie miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych wartości można odwołać się przez funkcję complete.cases().

# zastąp tylko te wartości, które nie są NA
vec[complete.cases(mat)] = mdl$cluster

W ostatnim kroku skopiujemy metadane obiektu landsat, ale tylko z jedną warstwą i przypiszemy mu wartości wektora vec.

clustering = rast(landsat, nlyrs = 1, vals = vec)

Sprawdźmy teraz jak wygląda klasteryzacja naszego obszaru na mapie.

colors = c("#086209", "#fdd327", "#d9d9d9", "#29a329", "#91632b")
category = c("lasy/woda", "pola uprawne", "odkryta gleba", "roślinność", "nieużytki")
plot(clustering, col = colors, type = "classes", levels = category)

Jeśli wynik jest satysfakcjonujący, to możemy go zapisać używając funkcji writeRaster() i wczytać w innym programie.

writeRaster(clustering, "clustering.tif")

Wizualizacje

Do analizy właściwości klastrów, zamiast statystyk opisowych, mogą zostać wykorzystane wizualizacje. Największe możliwości dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”.

ggplot2 wymaga przygotowania zbioru danych do odpowiedniej postaci. Dane muszą być przedstawione w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej zmiany można dokonać w prosty sposób używając pakietu tidyr.

# install.packages(c("tidyr", "ggplot2"))
library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych

Nasz zbiór danych jest całkiem pokaźny (blisko 9 mln wartości). Nie ma potrzeby przedstawiania wszystkich danych na wykresie. Wymaga to więcej RAM i znacząco wydłuża czas rysowania. Prawie identyczny efekt można uzyskać wykorzystując mniejszą próbkę danych. Jako przykład zobrazujmy jedynie 10 000 wartości z każdego kanału spektralnego. Do stworzenia losowej próby służy funkcja sample(). W wyniku otrzymamy indeksy wylosowanych wierszy.

idx = sample(1:nrow(mat_omit), size = 10000)
head(idx) # wyświetl 6 pierwszy indeksów
## [1] 294762 392686 640775 538191 270373 549593

Połączmy teraz wylosowane wiersze z macierzy z odpowiednimi klastrami (cbind()). Następnie macierz zamienimy na ramkę danych (as.data.frame()).

stats = cbind(mat_omit[idx, ], cluster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)
##          B1        B2        B3        B4        B5        B6        B7 cluster
## 1 0.0215950 0.0250875 0.0423575 0.0363350 0.1958350 0.1312650 0.0719200       1
## 2 0.0209625 0.0296525 0.0567400 0.0480775 0.3601475 0.1481500 0.0784650       4
## 3 0.0611400 0.0707100 0.1034350 0.1209800 0.3034425 0.3146350 0.2814150       3
## 4 0.0356200 0.0396350 0.0560250 0.0620475 0.1954225 0.2082925 0.1317050       1
## 5 0.0338875 0.0392775 0.0642200 0.0594900 0.3336650 0.2240500 0.1200450       4
## 6 0.0487650 0.0554750 0.0807475 0.0867425 0.2840550 0.2358200 0.1891525       5

Jak można zauważyć, powyższe dane mają formę szeroką (każdy kanał spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu wykorzystamy funkcję pivot_longer().

stats = pivot_longer(stats, cols = 1:7, names_to = "band", values_to = "value")

Dla formalności możemy jeszcze zmienić typ liczbowy danych (klastrów i kanałów) na kategoryczny.

stats$cluster = as.factor(stats$cluster)
stats$band = as.factor(stats$band)

Struktura danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.

ggplot(stats, aes(x = band, y = value, color = cluster)) +
  geom_boxplot()

Zmieńmy kilka domyślnych parametrów żeby poprawić stylistykę ryciny.

ggplot(stats, aes(x = band, y = value, color = cluster)) +
  geom_boxplot(show.legend = FALSE) +
  scale_colour_manual(values = colors) +
  facet_wrap(vars(cluster)) +
  xlab("Kanał") +
  ylab("Odbicie") +
  theme_light()

Zmieniając facet_wrap(vars(cluster)) na facet_wrap(vars(band)), zamiast zestawienia kanałów w poszczególnych panelach, możemy zestawić klastry.